// Sample.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include <stdexcept>
#include <iostream>
#include <map>
#include <ctime>
#include <cstdlib>
#include "string.h"
#include "blinalg.hpp"
#include "superelement.hpp"
#include "mesh.hpp"



int _tmain(int argc, _TCHAR* argv[])
{
	myindex N = 3;
	int max_refine = 2;
	std::cout << "Solve the 1D Elliptic problems: -(pu')'+qu=f by FEM" << std::endl;
	
	myindex nr_elements;
//	read and refine mesh
	nr_elements = readMeshFE("fem");
	for(int i=0;i<max_refine;i++)
		nr_elements *= 2;

//	initialize superelements
	std::vector<FElement> elem(nr_elements);
	InitialMeshFE(elem,"fem");
	RefineMeshFE (elem, max_refine);

	std::cout << "Compute Galerkin matrices for FEM" << std::endl;
	FullMatrix<double> MassMatrix(N,N), StiffnessMatrix(N,N);
	MassMatrix[0][0]=1.0;
	MassMatrix[1][1]=2.0;
	MassMatrix[2][2]=3.0;
	MassMatrix.print();
//	computeStiffnessMatrix (StiffnessMatrix, elem);

	std::cout << "Compute vector right hand side for FEM" << std::endl;	
	EVector<double> data(N), RHS(N);
	RHS.setZero();
	RHS[0] = 1.0;
	data[0] = 2.0;
	data[1] = 1.0;
	data[2] = 3.0;
//	data *= 2.0;

	data -= RHS;
	data.print();
	RHS = MassMatrix*data;
	RHS.print();
	std::cout << "Solve the linear system" << std::endl;	

	std::cout << "Compute the L2 error" << std::endl;	
	
	//Tinh ve phai

	return 0;
}

